Self-similar Bumps and Wiggles: Isolating the Evolution of the 
BAO Peak with Power-law Initial Conditions 
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Motivated by cosmological surveys that demand accurate theoretical modeling of the baryon 
acoustic oscillation (BAO) feature in galaxy clustering, we analyze N-body simulations in which 
a BAO-like gaussian bump modulates the linear theory correlation function ^z,(r) = (ro/r)"^"^ of 
an underlying self-similar model with initial power spectrum P{k) = Ak" . These simulations test 
physical and analytic descriptions of BAO evolution far beyond the range of most studies, since we 
consider a range of underlying power spectra (n = —0.5, —1, —1.5) and evolve simulations to large 
effective correlation amplitudes (equivalent to as =4—12 for rbao ~ 100/i^^Mpc). In all cases, 
non-linear evolution flattens and broadens the BAO bump in ^(r) while approximately preserving its 
area. This evolution resembles a "diffusion" process in which the bump width (Jbao is the quadrature 
sum of the linear theory width and a length proportional to the rms relative displacement Epair('"bao) 
of particle pairs separated by rbao- For n = —0.5 and n = —1, we find no detectable shift of the 
location of the BAO peaJi, but the peak in the n — —1.5 model shifts steadily to smaller scales, 
following rpeak/j'bao = 1 — 1 .08(ro /rbao ) ^' ^ . The perturbation theory scheme of McDonald (2007) [Jl 
and, to a lesser extent, standard 1-loop perturbation theory are fairly successful at explaining the 
non-linear evolution of the fourier power spectrum of our models. Analytic models also explain why 
the ^(r) peak shifts much more for n = —1.5 than for n > —1, though no ab initio model we have 
examined reproduces all of our numerical results. Simulations with Z/box = lOrbao and Lbox = 20rbao 
yield consistent results for ^(r) at the BAO scale, provided one corrects for the integral constraint 
imposed by the uniform density box. 
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I. INTRODUCTION 



The detection of the baryon acoustic oscillation (BAO) 
signature from observations of galaxy clustering 0, Q 
represents a triumph of large-scale-structure theory and 
of state-of-the-art cosmological surveys. The BAO fea- 
ture, imprinted by sound waves that propagate in the pre- 
recombination universe Q, provides a "standard ruler" 
that can be used to measure the distance-redshift rela- 
tion and the evolution of the Hubble parameter H{z) 0- 
0]. BAO measurements in the Sloan Digital Sky Survey 
(SDSS) yield a 2.7% measurement of the comoving dis- 
tance to z = 0.275 ([1,[1]; improved from the 4% precision 
of Several ongoing experiments - WiggleZ [lol [lli 
HETDEX [12, and the BOSS survey of SDSS-III [13 
- seek to extend these measurements to higher redshift 
and improve their precision, using spectroscopic surveys 
of galaxies and (in the case of BOSS) the Lya forest. 

seek 



Pan-STARRS [IJ and the Dark Energy Survey [If 
to measure the distance-redshift relation using the BAO 
feature in angular galaxy clustering, and the Large Syn- 
optic Survey Telescope [l6j will eventually reach much 
higher precision measurements. Other ambitious experi- 
ments - the ground-based BigBOSS survey and the 



space-based WFIRST [l8[ and Euclid [l9| missions - plan 
spectroscopic surveys of > 10* galaxies that in principle 
allow BAO measurements at the 0.1% level. 

The high anticipated precision of these experiments 
places stringent demands on theory. To fully exploit 
these measurements as probes of cosmic acceleration, one 
must understand the effects of non-linear gravitational 
evolution and non-linear bias of mass tracers (e.g. galax- 
ies or the Lya forest) on the location of the BAO feature, 
calculating any shifts to an accuracy below the statis- 
tical measurement errors. This challenge has inspired 
many analytic and numerical investigations of BAO evo- 
lution [20l - |29j , most of them focused on a ACDM cos- 
mological model (inflation and cold dark matter with a 
cosmological constant) with parameters close to those fa- 
vored by recent observations. In this paper, we pursue a 
complementary approach, inspired by N-body studies of 
self-similar cosmological models with a scale-free initial 
power spectrum P{k) = Ak^ [e.g. [30l - [37j . Specifically, 
we investigate models in which the correlation function 
of the initial density field (the Fourier transform of its 
power spectrum) is 
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a power-law modulated by a Gaussian bump centered 
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at a "BAO" scale rbao-^ For specified values of n and 
the bump height and width (Abao and tTbao)j the non- 
linear evolution of these initial conditions should depend 
only on the ratio tq /^bao of the correlation length to the 
BAO scale, and not (except for the overall change of 
scale) on the individual values of ro and rbao- Strictly 
speaking, this statement holds only for a particular cos- 
mological model (e.g. 57^ = 1, ^l\ — 0) in which the 
expansion factor a(t) is a powerlaw of time, but we will 
show that the bump evolution is nearly identical for an 
flm = 0.3, = 0.7 cosmology when evaluated as a func- 
tion of the linear growth factor. 

There are several valuable aspects of this approach. 
First, by varying n, abao and ro/rbao, we can investigate 
the interplay among power spectrum slope, bump width 
and non-linearity in determining the shape and location 
of the BAO feature. Second, we can test analytic (e.g. 
perturbation theory) descriptions of BAO evolution over 
a much wider range of conditions than they have been 
tested to date, to see how well they capture the underly- 
ing physics of BAO evolution as opposed to working in a 
specific case. Among other things, we evolve our simula- 
tions to values of tq /^bao much larger than those of con- 
ventional ACDM, so that we can clearly see where pcr- 
turbative approaches break down and how far they can 
be pushed. In this regard, our approach is similar to that 
of ^ and lU who use a "crazy" CDM (cCDM) model 
with parameters {flm = l,rif, = 0.4, erg = 1) designed 
to produce larger BAO wiggles and stronger non-linear 
effects. Third, the self-similarity of our model allows for 
numerical tests where, as a consistency check, the evolu- 
tion of the bump from simulations with the BAO bump 
with different numerical choices (e.g., box size relative 
to BAO scale, mean interparticle spacing, time steps) 
should all agree when compared at the same ro/rbao- 

Qualitatively, one expects the non-linear evolution of 
the BAO feature to involve a broadening and attenuation 
of the bump in configuration space, as discussed by j2^ . 
who describe matter scattering out of the BAO "shell" . 
In Fourier space this phenomenon is seen as a damping of 
oscillations at high-A:. In many pcrturbative approaches, 
this damping is exponential with a scale, S, given by 



1 



(2) 



For pure powerlaw cosmologies one can easily see that 
this expression will be problematic. Physically, Eq. [2] is 
the rms displacement of particles - which includes the 
contribution from bulk motions that shift all particles in 
a large volume coherently - whereas the damping of the 
BAO feature is more fundamentally related to the rms 
relative displacement of pairs of particles. For the mod- 
els investigated in this paper this subtlety becomes very 
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^ Note that a pure power-law spectrum P(k) 
to a correlation function ^(r) oc j-— [Jg 



Ak" corresponds 
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FIG. 1: A comparison of the linear theory matter autocor- 
relation function for ACDM (black, becoming dashed when 
< 0) and the linear theory matter autocorrelation func- 
tions investigated in this study. The ACDM correlation func- 
tion shown was generated using the fiducial WMAP7 cos- 
mology (assuming flatness), and the amplitude shown cor- 
responds to 2 = 0. For comparison these different clustering 
distributions are normalized to have the same non-linear scale, 
ro, as the ACDM case, where ^^(ro) = 1. 



important, and we argue that the broadening of the bump 
in our simulations scales according to the rms pairwise 
displacement equation (Eq. [TT] below) . 

We describe our initial conditions and simulation setup 
in § [ni show and characterize our results for the bump 
evolution in § HID and establish the numerical reliability 
of our results with self-similarity tests in § IIVI In §|V]we 
show the power spectra in our simulations and compare 
both phenomenological and ab initio quasi-linear models 
to the simulation results. We compare our results with 
this setup to canonical ACDM in § IVII and comment on 
the broader relevance of our findings. Finally in § IVIII 
we summarize our main conclusions and mention future 
directions for investigating this model. 



II. SIMULATIONS 
A. Initial Conditions 

We generate the initial conditions for the simulations 
by fourier transforming Eq. [1] to a power spectrum, 
Pic(fc), and using the pubhcly-availablc code 2LPT [39| . 
which computes particle displacements with second-order 
Lagrangian perturbation theory, to generate particle ini- 
tial conditions files. 2LPT has been shown to minimize 
transients compared to the first order [iO] approximation. 

In Fig. [1] we compare the three different ^ic('") mod- 
els explored in this paper (blue, green, and red) to a 
standard ACDM correlation function (black). We show 
the fourier transform of these correlation functions - 
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n = -0.5 
-n = -1.0 
-n = -1.5 
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TABLE I: Normalization/Conversion Table 
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FIG. 2: The linear theory power spectra of the models shown 
in Fig. [T] with the same normahzation. 



the resultant Pic(fc) - in Fig. [2] compared to a flat 
ACDM power spectrum generated from CAMB [4l| as- 
suming fiducial WMAP7 parameters [l^. In keeping 
with convention, we refer to the powerlaw in fourier 
space {n = —0.5, —1.0, —1.5) rather than in configuration 
space. These choices for the powerlaw slope are inspired 
by the resemblance to the ACDM correlation function 
on different scales. Similarly, unless otherwise noted, we 
choose fTbao = 0.075 rbao as the ACDM-inspired gaussian 
width and Abump ~ 2.75 as the gaussian amplitude of 
the BAO feature. 

In this study our time variable is ro/rbao, where 
^i(ro) = 1. This quantity grows as the amplitude of 
^i,(r) becomes larger and the correlation length rg in- 
creases. For convenience we show conversions between 
this convention for the non-linear scale and other choices 
in Table ID Other popular conventions define the non- 
linear scale as cr(i?*) = 1 or CT(i?*) = Sc, i.e. the scale 
where the rms density in spheres reaches one or reaches 
the threshold for spherical collapse, 5c — 1.69. We show 
R*/^ha.o for cr(i?,) = 1 in the second column in Ta- 
ble HI to convert from <7{R^,) = 1 to a{R^,) = 5c, mul- 
tiply this column by 5^^"^^^^ ■ A fouricr-space conven- 
tion for the non-linear wavenumber, A^(fcNL) = 1 where 
A^(fc) = k^P{k)/{2n)'^, is also shown in the third col- 
umn, /cnl is shown divided by fcbao ~ 27r/rbao so as to 
be independent of a specific choice of rbao and to reflect 
the self-similar nature of the setup. Finally, the fourth 
column shows the effective value of cts , computed assum- 
ing Tbao = 100/i~^Mpc. More generally this column can 
be interpreted to be the rms density contrast in spheres 
of radius 8% of the BAO scale. 

We begin our simulations at the earliest epoch listed 
for each of the three models shown in Table HI and we 
obtain outputs at each of the epochs listed. 
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B. Approximate Solution for Pic{k) 



Starting from the fourier transform relation. 



Pic(fc)=47r/ eicM ^^r^dr, (3) 
Jo kr 

and breaking up ^ic{r) in Eq.[T]into two terms, we expect 
Picik)^ PpoAk) + Pwie{k). (4) 
An exact analytic solution exists for the powerlaw term 



the fourier transform of P, 



pow 



(ro/r)""'"'^ with amplitudes related by 
27r2 {2 + n) 



Aa^k"" is ^(r) 



n+Z _ A n+3 
— ^ri ' 



r(3 n) sin((2 n)7r/2) " 
The remaining Pwig(fc) term in Eq. |3]is given by 



(5) 



^wig(^) ~ 



47rAbump?'o 



-(«+2)g-(r-rb„)-/2fT^,, 



sin(A:r) dr 



(6) 
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Up to a normalization, the integral is simply the expec- 
tation value of 7-^("+2) sin/cr over a gaussian probability 
distribution p{r) centered on rbao with width Ubao (but 
truncated at r > 0): 

r-("+2)sin(/cr)p(r) dr w (27rag^Ji/2( r-("+2) sin(fcr)). 

(7) 

Since p{r) is strongly peaked at r = rbao, and since 
sin(/cr) is generally much more sensitive than 7-"^"+^) to 
the value of r,^ we have, to good approximation, 

(r-("+2) sin(fcr)) w (r-("+2)) (sin(fcr)) « r~'-''+^\sm{kr)) , 

(8) 

leaving only the expectation value of sin(fcr) to be deter- 
mined. This expression is given by 

/"CxD 

( sin(fcr) ) = (27rcrg^J-i/2 / ^-{r-r^.a^/^-yLa sin(fcr) dr 

Jo 

w sin(fcrbao)exp(-(fcfTbao)^/2). (9) 
This line of approximation ultimately leads to 

Fic(fc) « A„r3(fcro)" + (10) 

25/2^3/2. 2 / rp \"+' sin(fcrbao) fc2,g^ 

^ ^' ^bump^bao' bao ' I , 



rba 



With our ACDM-inspired choices for the constants in 
this expression (discussed in ^ III A[) . our approximation 
for Pic (fc) agrees with the numerical integration to better 
than a percent (relative to the underlying powerlaw) over 
the entire range of fc- values. 



C. Integration of Particle Trajectories 

We used the publicly-available Gadget2 code [i^] to 
integrate particle trajectories from the initial conditions. 
Gadget2 is a hybrid, Trce-PM code in which the long- 
range gravitational forces are computed by solving the 
Poisson equation in fourier space while the short range 
forces are computed using a Tree algorithm [U- Gad- 
get2 is parallelized using standard MPI and allocates pro- 
cessors/cores with the space-filling Peano-Hilbert curve. 
This allows the code to perform well on massively-parallel 
machines. 

Throughout, unless otherwise noted, we simulate the 
powerlaw times a gaussian model using a flat r2,„ = 1.0 
cosmology with no dark energy, much like in self-similar 
pure powerlaw investigations [e.g. HO, H^] or in cCDM 



^ sin(fcr) goes as when k is small, and clearly varies rapidly 
with r when k is large. By contrast varies as r~^-^ for 

n = —1.5 and r~^-^ for n = —0.5. Most of the inaccuracy in the 
final result for Pic{^) comes from Eq. [5] The approximations in 
Eqs. [7] & E] are more accurate because they only depend on the 
assumption that cxp{—r^/2a^^^)dr^0. 



[20l l2l| . This choice allows structure to grow indefi- 
nitely, avoiding the freeze-out limit when the dark en- 
ergy component comes to dominate. However, in § IIVB| 
we present some simulations that include a cosmological 
constant and conclude that the evolution of the bump 
still only depends on the ratio of the non-linear scale to 
the BAO scale, even when dark energy is present. 

Most of the simulations presented here, unless other- 
wise noted, were run with 512'^ particles using a 768'^ 
PM grid for the large scale forces and a comoving force 
softening (relevant to the tree part of the code) of 1/ 4th 
the initial mean interparticle spacing. Our box size was 
chosen to be ~20x larger than the BAO scale, making 
the force softening ~l/2000th the scale of the box. We 
ran seven realizations of each model in order to obtain 
better statistics on large scales. We also performed pure 
powerlaw simulations (i.e. no wiggles) with the three 
cases (n = —0.5, n = —1, n = —1.5) to compare with the 
cases that include a BAO feature f Appendix E)) . Also 
note that we apply a correction to ^mcas('') to account 
for the artificial enforcement of the integral constraint 
on ^moas(f) (Appendix [B|) . This correction is important 
on large scales for n < —1. 

The simulations were evolved to the point where the 
non-linear scale reached approximately 30% of the initial 
BAO scale. As in pure powerlaw simulations, there is a 
concern that for steep power spectra the missing power 
on scales larger than the box will invalidate the results. 
However, even in the last output of the n = —1.5 case, 
which has the most large scale power, our simulations fall 
well within the guidelines recommended by [s^ , and the 
self-similarity of the pure powerlaw results in Appendix 1X1 
seem also to confirm the validity of our simulation results. 

All of the simulations presented here were performed 
using the Glenn cluster at the Ohio Supercomputer Cen- 
ter'^. In total, the results in this paper are based on 28 
512'^-particle simulations of powerlaw-|-bump initial con- 
ditions, 21 256'^-particlc simulations and 28 512'^-particle 
simulations used in the tests of IIVI and 20 512'^-particle 
simulations of pure powerlaw models presented in Ap- 
pendix 13 



III. EVOLUTION OF THE BAO BUMP 

A. ^(r) results for fiducial case 

Fig. [3] presents our main results for the configuration- 
space evolution of the BAO feature. Remarkably, when 
divided by the pure powerlaw correlation function as in 
the plots on the right hand column, the BAO feature 
maintains a gaussian shape throughout the non-linear 
broadening and damping that occurs in structure for- 



' |http : //www ■ osc ■ edu/supercomputlng/hardware/ 1 



5 




bao 





FIG. 3: Matter 2-point correlation function results for a powerlaw times a gaussian model for dark matter clustering. The 
upper two panels show results for the n = —0.5 background powerlaw, while the middle two panels show n = —1 and the 
lower two n — —1.5. The left panels show the measured matter autocorrelation function from the simulations at various epochs 
as colored points and, in dashed lines with the same color scheme, the corresponding linear theory correlation function at 
each epoch. The right panels show the matter autocorrelation function divided by the pure powerlaw correlation function, 
^pow(r). Black dashed lines show the linear theory prediction, which is independent of epoch. Typical errors on the mean 
for ^incas{r) /(,pow{r) are shown off to the right for various epochs, but note that errors are strongly correlated across the full 
range of the bump. On the right hand panels we also overplot with dot-dashed lines the best fit gaussians with the same color 
scheme as the measurements from simulations, ^(r) has been corrected for the integral constraint as described in Appendix IB] 
Comparable plots for a ACDM power spectrum appear in Fig. [15] 
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mation.^ In linear theory the bump would maintain the 
initial shape as indicated with the black dashed lines on 
the right hand column. 

We overplot the best fit gaussians on the right hand 
column with dot-dashed lines of various colors corre- 
sponding to different epochs to emphasize and better il- 
lustrate this gaussian behavior. We consider quantitative 
measures of the evolution in bump amplitude and width 
in§|llICl 

When comparing the three models, one should bear in 
mind that at fixed rp/rbao the bump in the n = —0.5 
case is at a much lower clustering amplitude than in the 
n = —1.5 case simply because an n = — 1.5 powcrlaw has 
much more large scale power, and we defined the initial 
bump feature to be a gaussian times (rather than added 
to) a powerlaw. The simulation data for the n = —0.5 
case are noisier, especially at early epochs, because we 
are measuring a weaker signal. 

The other striking feature of Fig. |3]is that the location 
of the bump maximum stays nearly fixed in the n ~ —0.5 
and n = — 1 cases, even when they are evolved to high 
values of ro/rbao (corresponding to cts = 6 — 12), while 
the location of the maximum for the 7i = — 1.5 case shifts 
substantially at late times. The shifts for n = —1.5 are 
6, 14, and 29 % at ro/rbao = 0.153, 0.263, 0.386 (corre- 
sponding to (Tg = 2, 3, 4). By contrast, in ACDM one 
typically sees shifts of ^ 0.5% by z = (cts ~ 0.8), and 
extrapolating the fitting formula of [l^l to an extreme 
value of (Tg — D{z)/ D{0) w 4 predicts a shift of only 
~ 5%. Qualitatively, we can understand the different be- 
havior of n = —1.5 as a consequence of the much higher 
clustering amplitude at r « rbao (see Fig. [1]). We will 
discuss the non-linear shift of the BAO peak in further 
detail in following sections. 

As one last qualitative note on the n = —1.5 results 
in Fig. [21 at the two latest epochs one can see that the 
correlation function at r ^ 0.5rbao is showing significant 
non-linear evolution away from the initial power-law, in 
contrast to the other two cases. We avoid this region in 
determining the best fit gaussians to the simulation data. 



B. Evolution of a "Skinny" Bump 

We also investigated a case where the initial gaussian 
width of the bump was half of the value in the fidu- 
cial case, i.e. dbao = 0.0375 rbao instead of the ACDM- 
inspired value of dbao = 0.075 rbao- Keeping Abump fixed 
at 2.75, we performed simulations only for the n ~ —1 
background powerlaw. These results are shown in Fig. |4l 
The bump clearly maintains a gaussian shape as it is 
damped out, and, as in the fiducial n = —1 case, there 



* The exception, discussed below, is at late times (high clustering 
amplitudes) in the n = —1.5 model. 



does not seem to be any shift in the BAO peak by the 
end of the simulation. 



C. Quantitative Characterization of the Bump 
Evolution 

In Figs. [5] and |6] we plot evolution of the amplitude, 
width, and peak location measured from gaussian fits to 
our simulation results. These gaussian functions were 
determined by first making a rough determination of 
the BAO peak and bump amplitude from C(r)/^pow(r), 
then varying rbao, ^bump and tTbao in a 3-dimensional 
to find the best fit. This minimization was done using 
^(r)/^pow(r) as in the right hand panels of Fig. |3] rather 
than ^(r) itself. We avoided correlation function data 
more than Ar ^ 0.3rbao below the peak in finding the 
best fit gaussian, to avoid effects of non-linear evolution 
of the underlying powerlaw correlation function. Error 
bars in Figs. [S] and [5] were determined via jackknife er- 
ror estimation by sequentially omitting the correlation 
function results for one of the seven realizations and de- 
termining the best fit gaussians in each case. The errors 
on ^bump X Cbao/'^bao, a dimcnsiouless proxy for the area 
of the bump, are from propagated errors in the values of 
^bump and (Tbao- The n = —1.5 case suffers from a slight 
degeneracy between the amplitude of the bump and the 
magnitude of the non-linear shift, so the jackknife error 
bars arc slightly larger in this case. 

Fig. [5] shows our main results for the quantitative evo- 
lution of the dimensionless bump width, dbao/rbao, bump 
amplitude, Abump, and area Abump x (Tbao /rbao- In the 
top panel, in all cases there is significant broadening of 
the bump, while in the middle panel, even apart from the 
dot-dashed models which will be discussed in a moment, 
the amplitude of the n ~ —0.5 case appears to decrease 
more slowly than that of the other setups. 

The lower panel of Fig. [5] shows that the area under 
the bump stays remarkably constant, closely following 
the black horizontal dashed and dot-dashed lines as the 
bump broadens and attenuates. We speculate that the 
non- linear dynamics of the growth of structure is just dif- 
fusively moving apart the pairs at separation r ^ rbao so 
that (T^g^Q « (T^Q + (T^jg, where die is the initial bump 
width and a^is is the rms broadening from this diffusion 
process, while the area under the bump stays constant 
and the gaussian shape is maintained. These assump- 
tions underlie the models plotted in the top two panels 
of Fig. [5l The broadening is modeled by identifying a^^g 
with the linear theory equation for the mean-squared rel- 
ative displacement between pairs [Eq. 9 from |22| . 

f°° k^Hk 

^U^rlJ P{k)f^^{kn,), (11) 

Jo ^T'' 

where ri2 is the separation and 

2/1 sin(a;) 2cos(x) 2sin(a:) ^^ 
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FIG. 4: Results for a setup where the initial gaussian width of the bump is half of its fiducial value and the background 
powerlaw is set to n = — 1. As in Fig. [3] different outputs are shown in dilTerent colors, with the typical errors on the mean 
offset to the right, ^(r) has been corrected for the integral constraint as described in Appendix iBl 



In the limit ri2 — > oo, Eq. [TT] reduces to Eq. [21 i.e., 
the rms pairwise displacement Epair asymptotes to the 
Zel'dovich displacement. However, modes with fcri2 1 
move pairs of particles separated by ri2 coherently, and 
while these modes may dominate the "bulk flow" they 
cannot affect clustering on scales < ri2 . Notably, Eq. [5] 
is infrared divergent for n < — 1 , while Eq. [TT] is IR con- 
vergent for n > — 3, failing only when the density contrast 
(not peculiar velocity) has a divergent large scale contri- 
bution. 

If, as an approximation to our model, we consider a 
pure powerlaw power spectrum, P{k) « Aa^k^, Eq. [TT] 
can be re-written as 

Spair=^^/ ^'^V\\i^)dx, (13) 

^~ Jo 

where x = kr 12- Selecting ri2 = ^bao, and utilizing 
Aa?' ~ fo^"^ (Eq. [5]), this implies a scaling of the form 

^.— rlJ^y. (14) 

\ ' bao / 

The width of the bump can therefore be modeled with 

2 2 I Y^2 

"ha.o — "IC + ^pair' 

2 _ 2 2 ( TO (15) 

\ 'bao / 

We use the symbol k„ and include a factor of 2 to em- 
phasize our characterization of the bump evolution as a 
diffusion process. For — 3 < n < —1, evaluating the inte- 
gral in Eq. [T3] yields an ah initio prediction for k„ (from 

S^ai.) Of 

_ 2 + nT{\ + n) sin(?i7r/2) 

~ 2-nr(3-|-n)sin((2 + n)7r/2)' ^ ' 



For shallower power spectra, n > — 1, both Eq. [13] and 
Eq.[2]are UV divergent and undefined. In Fig. [S] we 
therefore model the evolution of Ubao for n = —0.5 and 
ri = — 1 by assuming the scaling in Eq. 1 141 and empirically 
fitting K„ to our simulation results. For n ~ —1.5, Eq. 1161 
yields k_i.5 = 4/7. 

The curves in the upper panel of Fig. [S] compare Eg. 1151 
with values of k„ = {4.5,1.3,4/7} for n = —0.5,-1 
and —1.5, respectively, to the measurements from simu- 
lations. The Kn values for n = —0.5 and —1 were chosen 
by a visual fit to the simulation points, while K-1.5 is 
predicted ab initio. The model provides a good match to 
the data for ro/rbao < 0.1. Most significantly, the same 
Kn fits both the fiducial and skinny n = — 1 cases, sup- 
porting the conjecture that the bump width is effectively 
set by a quadrature sum of the linear theory "intrinsic" 
width and the rms pairwise displacement (Eq. I15p . The 
scaling with rms displacement holds fairly accurately out 
to large rp/rbao- For comparison with an analogous test 
in ACDM, Fig. 3 of [H] shows that Eq. [TT] accurately 
predicts the rms displacement of pairs initially separated 
by r ~ 100h~^ Mpc. This agreement extends to late 
times where the rms displacement has reached ^ 8% of 
this initial separation. The n = —1.5 results in the up- 
per panel of Fig. [5] suggest that Eq. [TT] is accurate (i.e. 
within errors) for rms displacements as large as ^ 15% 
of the scale of the initial separation. 

The bottom panel of Fig. [5] shows that the constant- 
area approximation holds well for n ~ —I and n ~ —1.5, 
but it breaks down for n = —0.5 when ro/rbao ^ 0.1. 
This lack of constant-area behavior for n = —0.5 at late 
times explains the divergence of points and model curve 
in the middle panel of Fig. [S] 
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FIG. 5: Results for the dimensionless width (top), amphtude 
(middle), and a proxy of the dimensionless area under the 
bump (lower panel). Overplotted in the top two panels is a 
diffusion model in which the broadening of the width scales as 
suggested by the rms pairwise displacement equation (Eq. Ilip 
while the area of the gaussian bump is held constant. The 
diffusion constant is predicted ab initio for n — —1.5 and 
chosen by visual fit to the data points for n — —1 and n — 
—0.5. Error bars are from jackknife error estimation. 



D. Movement of the BAO peak 

Fig. ini shows the change in position of the bump maxi- 
mum, determined as described in ^ IIII Cl bv fitting a gaus- 
sian to the ratio of the non-linear correlation function to 
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FIG. 6: Results for the non-linear shift of the BAO peak 
measured from the right panels of either Fig.|3]or Fig. 3] The 
"skinny" bump results are offset to the left by Aro/rbao = 
0.01 so as to avoid overlap with the fiducial n = —1 results. 
Dashed lines show a prediction for the shift of the peak based 
on Eq. 32 from [2^. Error bars are from jackknife error esti- 
mation. 



£ (r) = A_dUr)/dr^(^>(r) 



n = -1.5 




FIG. 7: Residuals showing the non-linear shift from subtract- 
ing the gaussian fits centered on the un-shifted BAO scale 
from the matter correlation function results from the first two 
outputs in the n — —1.5 case (blue and green points). The 
solid lines show the predictions of the [23] ansatz in this case 
using Amc ~ 34/21 for both outputs. The thin green-dashed 
line uses this ansatz but assumes a broadened and attenuated 
bump in Cmc(f) as discussed in the text. 



the linear-theory powerlaw. As already noted in our dis- 
cussion of Figs.|3]and|ll there is no significant shift of the 
peak location on our simulations for either the n = —0.5 
or n = — 1 cases (fiducial or "skinny" bump). Error bars 
on the n = —0.5 peak location become large at late times 
because the bump itself flattens and the large scale cor- 
relation is weak. For n = — 1, the skinny bump errors are 
initially lower than those of the fiducial model because 
the sharper peak can be centroided more precisely, but 
they are higher at late times because the skinny bump 
gets depressed to a lower amplitude. In contrast to the 
other cases, the n = —1.5 model shows significant and 
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strong peak shifts, evident already for rg/rbao = 0.024 
(cts = 0.5). Indeed, we have truncated the pfot before 
the final n = —1.5 output, with ro/j'bao = 0.386 and 

?'poak/?'bao = 0.71. 

We compare these results to an elegant model for the 
shift from [2^ that uses linear theory velocities and the 
pair-conservation constraint on ^(r) to track the average 
motion of pairs separated by rbao- Their equation (32) 
can be written 



Dl 



1 = 



3 dr 



(17) 



where D{z) is the linear growth factor, the subscript ic 
refers to initial conditions when fluctuations are fully in 
the linear regime, rbao is the linear theory BAO position, 
''peak is the non- linear position of the peak, and ^(r) is the 
volume-averaged correlation function interior to radius 
r. For D{z)/Dic ^ 1 and our initial conditions, this 
equation leads to the approximate result 



'"peak 
'"bao 



n+3l l/("+3) 



(18) 



where C„ would be 1.0 for a pure powerlaw spectrum and 
incorporating the bump gives C„ ~ {1.13, 1.26, 1.38} for 
n = {— 0.5, — 1, — 1.5}. For n < this formula predicts 
that the peak shifts to smaller scales. In the limit of 
small ro/rbao, a binomial expansion yields 



'"peak 
^bao 



n 



^bao 



n+3 



(19) 



Since ''g'*''^ oc D'^{z), the non- linear shift grows as the 
square of the linear growth function as expected from 
PT [e.g. [2l|, [l^l, and the displacement is larger for more 
negative n. 

The seemingly quite different argument of leads 
to a similar expression for the peak shift. They propose 
modeling the non-linear correlation function in the neigh- 
borhood of the bump by 



eNL(r)«e(r)+^ 



dr 3 



(20) 



where the mode-coupling factor can be treated as 
a fitting parameter but the value 34/21 obtained from 
PT is in fact close to the best-fit numerical value (see 
psf . Appendix A). With judicious use of Taylor expan- 
sions in the limit of small shift and minimal non-linear 
broadening, one can derive 



^pcak _ 2. -I- — — '"0 

^bao 

21 n Vrbao 



ri-l-3 



(21) 



hence a shift about 50% larger than Eq. ([T9| but with 
the same dependence on rg and n. 

Dashed lines in Figure [5] show the prediction of 
Eq. The model correctly predicts that the shift is 



much larger for n = —1.5 than for n = —1 or n = —0.5. 
For n = —1.5, it tracks the numerically measured shift 
remarkably well. For the other n values, it predicts too 
large a shift for rg > O.lrbao! at smaller rg, the model 
is consistent with the numerical results within the error 
bars, but the numerical results are also consistent with 
zero shift. We note that our treatment docs not include 



the 1-loop PT extension of [25|'s model, which could im- 
prove agreement at later epochs. 

Figure [7] compares Eq. ((20|) to the first two outputs 
of the n = —1.5 simulations. For rg/rbao ~ 0.024, this 
model predicts the distortion in the neighborhood of the 
peak remarkably well, with no free parameters. Note that 
there is a clear non-linear shift of the peak at this output, 
despite the low value of as = 0.5. For rg/rbao = 0.061, 
the model predicts too large a distortion. However, if we 
insert the broadened and lower amplitude bump (taking 
CTbao and Abump from the model discussed in the previous 
section) into the calculation of Eq. (PH]) , an approach that 
seems reasonable but is not rigorously justified, then we 
get the dashed green lines in Figure [3 which agrees much 
better (though not perfectly) with the numerical results. 

We conclude that these analytic approaches can ex- 
plain why the shift in the bump location is much larger 
for 71 = —1.5 and can capture at least some of the quan- 
titative behavior of the peak shift. However, they do 
not work accurately over a wide range of rg/rbao and n. 
We will return to the comparison of PT predictions and 
our numerical results in i fVl in the context of the power 
spectrum. 



IV. SELF-SIMILAR TESTS 

In an fi™ = 1 pure powerlaw model, i.e. P{k) = Ak"^, 
since the only scale germane to the problem is the am- 
plitude, A, the evolution of clustering statistics should 
depend only on the value of A or some derived variable 
such as fcNL = (27rV^)^/^"+^^- The evolution may be dif- 
ferent for each powerlaw but with n fixed there should be 
a unique function (e.g. of fc/fcNL, or r/_RNL, or M/Mj^i^, 
...) that fully describes any given clustering statistic, 
even well into the non-linear regime. In the early days 
of cosmological N-body investigations, demonstrations of 
self-similar evolution with pure powerlaw cosmologies, 
in addition to providing physical insight, also gave de- 
cisive confirmations of the accuracy of simulations [e.g. 
I30l.[3ll [33l - l35| . We take advantage of the simplicity of the 
powerlaw times a gaussian setup to perform self-similar 
tests that can be used in an analogous way to test the 
accuracy of the simulations on the scale of the bump. 

The powerlaw times a bump setup clearly has two 
scales at play instead of one, so in this case, for a given 
powerlaw and a given initial bump width, the non-linear 
dynamics should evolve only as a function of the ratio 
of the non-linear scale to the BAO scale. The dynam- 
ics are self-similar in the sense that any property of the 
system, such as the broadening of the bump or the shift 
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FIG. 8: Tests of robustness to numerical parameters. {Top ) Comparison of the bump region of ^(r) in simulations with 
'"bao/ibox = 1/10 (points) to the gaussian fits (lines) from the fiducial simulations in Fig. O which have rbao/^box = 1/20. 
{Bottom) Comparison of simulations with rbao/^box ~ 1/20 but 256'^ particles, hence rbao/^p ^'''^ ~ 12.5, to the fiducial 
simulations with 512"^ particles and rhuo/np ^^"^ — 25. 



in the peak for a particular powerlaw, is determined by 
how close the non- linear scale, tq, has come to the BAO 
scale. Unlike a ACDM simulation the result should not, 
in principle, depend on whether the bump is initially set 
at, e.g., 100 Mpc or 130 Mpc; only the ratio of 
the non-linear scale to the BAO scale matters in deter- 
mining the evolution. If the N-body results do depend, 
separately, on the BAO scale or the non-linear scale, this 
can be interpreted as a sign of numerical artifacts. 



A. Robustness to Varying Box Size and Mean 
Interparticle Spacing 

Cosmological N-body simulations unavoidably intro- 
duce two artificial scales into the problem - the box 
size, Lbox, and the initial mean interparticle spacing, 
Ip = rip ^^■^ = Lbox/A^^^"^- Both of these scales can po- 
tentially interfere with the evolution of the BAO feature 
and bias one's results. In the upper panels of Fig. [S] 
we show results from tests where the BAO scale has 
been doubled (or equivalently the box-size halved), such 
that Tbao/ibox ~ 1/10 instead of the fiducial value of 
fhno/Lhox ~ 1/20 in the simulations shown elsewhere 
in the paper. The number of particles in this test 
is kept fixed at 512'^, so that the ratio of the BAO 
scale to the mean interparticle spacing increases from 

— 1/3 

Thao/np — 25 (as in the fiducial simulations) to 50. 



We also show tests (lower three panels) where the box 
size is kept fixed while the number of particles is de- 
creased to 256'^, arguably more akin to a conventional 
convergence test. Each panel in Fig. [5] shows the results 
from seven realizations (as in the fiducial simulation set), 
and in each panel we plot the best-fit gaussians from our 
fiducial set of simulations (dot-dashed lines). Note that 
for the "double-the-bump" tests in the upper panels of 
Fig. ISl these simulations had to be run for much longer 
than in the fiducial case in order for the non-linear scale 
to approach the BAO scale, which had been placed at 
twice the fiducial separation. 

To the extent that the simulations in Fig. [5] match 
the fit from the fiducial set of simulations, the evolu- 
tion can be said to be self-similar and unaffected by the 
artificially-introduced numerical scales. For the double- 
the-bump tests, the results seem to match the fiducial 
simulations well. In this case, especially for n = —1.5, 
the integral-constraint correction to ^ (r) discussed in Ap- 
pendix [B] is critical. We interpret this agreement as an 
indication that rbao/^box ^ 1/10 is acceptable if one 
includes integral-constraint corrections. Note that the 
measured errors on the mean are larger for these tests, 
which measure the correlation function on scales closer 
to the box scale than in the fiducial simulations. These 
larger errors are consistent with expectations from Gaus- 
sian statistics in a finite volume j46| . 

The 256"^ test was not quite as successful. The acccl- 
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FIG. 9: Results for a model including a cosmological constant 
(yim = 0.3, Q.A = 0.7 at the output with ro/rbao = 0.043) and 
with an n = — 1 background powerlaw. The first and third 
outputs (blue and red) are directly comparable to the first 
and second outputs of the fiducial n = — 1 simulations; thus 
the gaussian fits to those outputs in the fiducial case are over- 
plotted. The second and last outputs (green and cyan) are 
compared to extrapolations from the fiducial n = — 1 case 
assuming no non-linear shift and a model for the bump evo- 
lution as described in § HIT CI 



erated attenuation of the bump in the n = —1.5 case 
is severe enough to be of particular concern, especially 
since this setup is the one which actually sees an ap- 
preciable change in the BAO peak. The n = —\ simu- 
lations agree much better but still shghtly underpredict 
the bump height. This also seems to be the case with 
the n = —0.5 results, which are more noisy. Though not 
quite a failure, we interpret this test to recommend keep- 



ing rbao/n 



-1/3 



> 25. as in our fiducial set of simulations. 



The tests in Fig.[5]show that the evolution of the bump 
- its flattening, its movement in the n = —1.5 case, the 
lack of movement in the n = —0.5 and —1 cases, and the 
unexplained behavior of the bump area in the n = —0.5 
case - is robustly predicted even when numerical param- 
eters are changed substantially. For the wider impor- 
tance of using BAO to constrain cosmology, this is an 
encouraging sign that modest N-body simulations can 
accurately render the non-linear shift of the BAO peak 
with very different models for the broad-band clustering. 
For power spectra that span a much wider range than 
ACDM models, numerical parameters 7'bao/-^box ^1/10 
and rbao/?^p ^ 25 appear to be adequate. 



B. A Test with Dark Energy 

Even with a powerlaw initial spectrum, the introduc- 
tion of dark energy in principle breaks self-similarity by 
defining a characteristic time (when Vim and Q\ are 
equal), however in linear perturbation theory and the 
quasi-linear Zel'dovich and adhesion [U I3 approxima- 
tions, evolution is determined only by the linear growth 



factor, with no direct dependence on Pm{o) or pDE(a)- 
[iot and [l^l demonstrate that this dependence on the 
Hnear growth factor alone remains a very good approxi- 
mation in fully non-linear N-body simulations, the latter 
also showing explicitly that the full equations of motion 
for cosmological perturbations are weakly dependent on 
the individual values of ilm and JIa when those equations 
arc expressed using the linear growth factor as the time 
variable. We may therefore expect self-similar evolution 
of the BAO bump as a function of ro/rbao, even when 
dark energy is included. 

Fig. [5] compares these expectations to the N-body sim- 
ulation results by presenting the evolved bump in a set 
of simulations (n = — 1) that include a cosmological con- 
stant in comparison to the gaussian fits to our fiducial 
simulations. The simulation set has VLm = 0.3, f^A = 0.7 
at output ro/rbao ~ 0.043. For some outputs, we have 
interpolated between outputs of our fiducial simulations 
assuming the model for the bump evolution discussed 
in § IIII CI A substantive difference between Q.^ = 1 
and dark energy models is that the growth of structure 
"freezes out" as dark energy becomes the dominant com- 
ponent of the universe. The last output in Fig. [S] is very 
close to this "freeze-out" limit in the linear theory growth 
function, which prevents ro/rbao from growing beyond 
0.073 in this case. 

The gaussian fits to the fiducial simulations agree well 
with the simulation results in Fig. |9l even for the last 
output which, with considerable computational expense, 
was evolved very close to the freeze-out limit. This con- 
firms the expectations of self-similar evolution for this 
setup even in cosmologies with dark energy. 



V. EVOLUTION OF THE BAO FEATURE IN 
FOURIER SPACE 

A. Power Spectrum Estimation 

Power spectra were determined for the n = —0.5,-1 
and —1.5 models by mapping the particles onto a 1024"^ 
grid using the cloud in cell (CIC) assignment scheme. 
Performing a discrete fast fourier transform on this grid 
yields 5{k) and fourier amplitudes P{k) ~ \5{k)\'^. The 
artificial smoothing introduced by the griding scheme is 
corrected for by dividing P(k) by the appropriate assign- 
ment function for CIC [5l[, and the corrected P(k) is 
binned in k to yield P{k). Following [s^ we do not in- 
clude any kind of shot noise correction [e.g [H, [U, and 
we follow their advice in trusting the computed power 
spectra only up to half the particle nyquist wavenumber, 
as indicated with black dotted vertical lines in Fig. [Till 
which presents our primary power spectrum results. The 
power spectrum up to this fc-value should be negligibly 
affected by the aliasing of the 1024'^ grid. Notwithstand- 
ing our conservative decisions in measuring P{k), we will 
argue in the next section that a simple phenomcnological 
model that draws on results from pure powerlaw Simula- 
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FIG. 10: Left column: Measured power spectra for the fiducial n — —0.5 (top), n = — 1 (middle) and n = —1.5 (bottom) 
simulations. The x-axis is shown normalized to the scale of the BAO feature, fcbao ~ 27r/rbao and the y-axis is likewise shown as 
a dimensionless quantity, P{k)/r't,a,o- There is no correction for shot noise; the shot noise level is indicated with dot-dashed lines. 
Right column: Results from dividing by the linear theory pure powerlaw. In both columns a phenomenological model (Eq. 1221 
solid colored lines) is compared to the simulation results. The scale corresponding to half the particle nyquist wavenumber is 
indicated with a vertical black dotted line. 
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tions (Appendix |X| allows our predictions to be extended 
to much higher k for the early outputs. 

We report power spectra throughout, normalizing the 
wavenumbers by fcbao = 27r/rbao, and giving the power 
spectrum amplitudes in terms of P{k) /r^^^- This re- 
flects the self-similar nature of the problem and allows 
more straightforward identification of the fc-values of var- 
ious nodes and anti-nodes. For technical reasons we 
throw out the measurements of the spectral power for 
k K. 27r/Lbox, which should be computed separately from 
measurements at higher k because of the different statis- 
tics of mode-counting near the scale of the box. The 
power on these scales is also inevitably noisy because of 
the small number of modes. 



B. Interpretation 

Ignoring the wiggles in Fig. [TU] for a moment and fo- 
cusing on the evolution of the overall shape of the power 
spectrum, the results are bracketed by the n = —0.5 
spectrum, which trails behind the linear theory power- 
law at high fc, and the n = —1.5 spectrum, which clearly 
outpaces the linear theory clustering prediction. This 
behavior is expected from perturbation theory [ssl Isg} , 
and the trend is more clearly shown in the pure pow- 
erlaw plots in Appendix |^ Physically, the behavior of 
the n = —0.5 power law is sometimes described as "pre- 
virialization" |57| . where on small scales the clustering 
power is so high that as halos form they pull away from 
the expansion of the universe and the non-linear power 
spectrum falls behind the linear theory prediction. For 
much steeper powerlaws, like n = —1.5 or the high k 
spectrum of ACDM, the trend is the opposite; clustering 
power is "transferred" from large scales to smaller scales. 
The n = — 1 spectrum lies between these two extremes, 
and its spectrum is above and below the linear theory 
prediction in different ranges (Appendix [X| . 

With this in mind we modeled the power spectrum re- 
sults with a phenomenological approach, treating sepa- 
rately the non-linear evolution of the pure powerlaw spec- 
trum and modeling the wiggles by coupling the analytic 
solution in Eq. [TU]with the diffusion model introduced in 
§ ImCl Thus, the model is 



i^phon(fc) Anrl{kroTfn{k/k^i^) + 



(22) 



25/2^3/2^/ 



bump '-''bao ^bao 



2 ( rp \"+^ sin(fcr;^J .2^,2^^^ 



where cr(^ao from Eg. 1151 and, as in IIII C[ the area un- 
der the bump is assumed to be constant, AJj^^^p a'^^^ = 
^bumpCic. For the n = —0.5 and n = —I cases we 
assume no shift of the BAO scale, r[^^^ = r^ao, while 
for n = —1.5 we set the BAO scale using r^jg^^/rbao = 
1 — 1.08(ro/r-bao)^'^, which is a good description of the 
motion of the peak in Fig. [51 For the pure powerlaw evo- 
lution we use non-linear fitting functions to pure power- 
law simulations, fn{k/k^-i), which are described in Ap- 



pendix [Xj In Fig. [To] we show the predictions of the phe- 
nomenological model with solid lines up to the fc-values 
where the fitting function is well determined by the pure 
powerlaw simulations. 

This model works surprisingly well in the n = —0.5 
case, given that the constant area approximation seems 
to break down in the later outputs (Fig. [5]). The first few 
outputs of the n = —1 and n = —1.5 cases are also well 
matched by Eq. [211 For these first few outputs the phe- 
nomenological models may actually be more trustworthy 
than the simulation measurements: at high fc the pure 
powerlaw spectrum dominates, and the non-linear fitting 
functions in this regime are defined preferentially from 
later outputs in the pure powerlaw simulations, which 
should be unaffected by transients from initial conditions 
or shot noise. 

If the phenomenological model can be trusted at high 
k, our results for the first output shown in Fig. [TUj can 
be extended to fc/fcbao ^ 30 for n — -0.5, k/k\,s_o ^ 600 
for n = —1, and k/k\jiio ^ 50 for n = —1.5. Assuming 
again that simulations can be trusted up to half the par- 
ticle nyquist wavenumber, this is analogous to running 
simulations for this setup with ~ 2400^, ^ 48000^ and 
~ 4000"^ particles respectively^, assuming the same box 
size as the 512^ simulations presented here. 

At small fc and late times there are significant devia- 
tions, however, between the phenomenological model and 
the n — —1 and n ~ —1.5 results. Those outputs have 
features, especially around k ^ fcbao, that seem to be un- 
accounted for in Eq. 1221 In the next section we compare 
the simulation results to expectations from perturbation 
theory. 



C. Comparison with PT predictions 

Because of the IR divergence of / P{q) dq for steep 
power spectra, perturbation theory schemes that use this 
term to renormalize the higher order expansions will ei- 
ther be intrinsically problematic for these setups or in- 
volve non-trivial cancellations that make the numerical 
evaluation of perturbation theory predictions much more 
difficuh. Standard 1-loop PT (a^.a. SPT) is still for- 
mally well defined for n > —3 |3^ . [ssj and so, like [s^] 
who explored pure powerlaw spectra, we show predic- 
tions for this approach and for the closely related "cou- 
pling strength" RGPT^ scheme from [l|. In principle, 
with sufficient care to deal with infrared divergences, the 
predictions from a number of other perturbation the- 



-1 case comes from the 



^ The extraordinary value for the n 
high k fitting function from [37{ . 

^ This name reflects the approach of this scheme in which the be- 
yond Unear-order terms are introduced with a coupling strength 
parameter, A, and the solution is evolved from linear theory 
(A = 0) to the full non-linear prediction (A = 1) (P. McDon- 
ald, private communication). 
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ory schemes could be compared with the simulation re- 
sults presented in this paper. We explore the predic- 
tions of SPT and "coupling strength" RGPT as two rep- 
resentative schemes that are reasonably-well studied [e.g. 
lli,i3Z,^. 

Although, the n = —0.5 and n = — 1 cases are IR con- 
vergent for J P{q) dq, they arc still UV divergent. As a 
matter of principle, we regard UV divergences as less seri- 
ous than IR divergences - it is no surprise that perturba- 
tive calculations break down on small scales where fluctu- 
ations are large. However, as a practical matter they are 
still problematic. Through separating the powerlaw and 
wiggle terms we can avoid some of the cutoff dependence 
of the SPT predictions; the predictions shown for the n = 
— 1.5 case should be completely cutoff independent, while 
n ~ —0.5 and n = — 1 results depend on the UV cut- 
off. In what follows we choose fc„iax/fcbao ~ 160, but our 
qualitative conclusions would be unchanged even if this 
high-fc cutoff were increased by a factor of two. (More 
precisely, if the cutoff were increased by a factor of two 
the comparison with simulation results in Fig. [TU] would 
be quite similar, and the evidence in Fig. [T^] that SPT 
overpredicts the damping in the n ~ —0.5 case would be 
stronger. Fig. [T3l would be unchanged.) To capture the 
physics of the problem fcmax should be significantly larger 
than the wavenumbers relevant to the BAO feature, i.e., 
fcmax ^ fcbao, 27r/(Tbao- The latter constraint is more im- 
portant, suggesting fcmax > 27r/CTbao ~ 13.3fcbao- [131 

choose fcinax ^Ny = T^N^^^^ / L^oyi to make their PT pre- 
dictions for powerlaw initial power spectra, but we feel 
that setting this upper limit according to the parameters 
of a simulation is a questionable thing to do when mak- 
ing ab initio predictions of non-linear behavior. Our PT 
predictions were calculated by modifying the publicly- 
available copter code from [201 better accommodate 
powerlaw cosmologies and our setup. 

The primary PT results are presented in Fig. 1111 At 
each output predictions are shown up to fcNL, roughly the 
scales where these schemes are expected to break down. 
Generally, "coupling strength" RGPT and SPT/SPT+ 
give fair-to-good predictions for the non-linear damping 
of the wiggles (as discussed below, SPT+ uses the non- 
linear fitting functions in Appendix |X] for the powerlaw 
evolution). The good comparison with the simulations 
for the n = —1.5 case suggests that the non- linear shift 
can be adequately captured by PT. An exception to this 
is clearly the SPT+ predictions for the n = —0.5 case, 
which seem to significantly overpredict the damping of 



the BAO feature. We discuss the SPT/SPT-f- predictions 
in more detail in the next section, breaking up the cal- 
culation into different "interaction" terms in an effort to 
gain insight into the non-linear physics. Predictions from 
"coupling strength" RGPT in Fig. [Til were not calculated 
by breaking up Pic{k) in this way, since this scheme does 
a much better job than SPT in predicting the evolution 
of pure powerlaw spectra 



37| 



D. SPT and SPT+ 

The 1-loop correction to the linear theory power spec- 
trum is given by [32| 



P{k) ^ PUk) + P22{k) + Pisik) 



(23) 



where 



P22{k) = 



and 
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fc3 



98 (27r) 
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- 2rx) 
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3r + 7x- lOrx^ 
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2rxY 



(24) 
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252 {27rf 
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drPdkr) — - 158 
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lOOr^ 



1 



(25) 

Notice that in P22{k) and Pi3{k) the linear power spec- 
trum appears twice, and as a result these terms increase 
in amplitude as the linear growth function to the fourth 
power. 

For pure powerlaw spectra, by including UV and IR 
cutoffs and using sufficient care to avoid the singularity 
in the denominator of the kernel in ^22(^)1 these integrals 
can be computed analytically [12, HI]- In principle, it 
may also be possible to obtain an exact solution for 1- 
loop corrections to the analytic expression for Picik) in 
Eq.[Tni but the complexity of the P22{k) kernel is difficult 
to overcome or approximate. 

To organize the calculation and for the most clarity in 
physical interpretation, we calculate the 1-loop correc- 
tions by treating separately the "interaction" ^ terms that 



arise from inserting Phik) 
in P22{k) and Pi3(fc), 



PpoAk) + P«.s{k) (Eq.lH) 
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FIG. 11: A comparison of power spectrum results from Fig. [10] with quasi- linear predictions from standard perturbation 
theory (SPT/SPT-I-, dashed hues) and the "coupling strength" RGPT scheme (dot-dashed lines) from SPT-|- treats the 
pure powerlaw evolution differently than SPT, using a fit to pure powerlaw simulation results (Appendix instead of the 
SPT prediction for the pure powerlaw evolution [s^. In each plot the a::-axis limits are set to include the low k measurements 
from simulations up to half the particle nyquist wavenumber, approximately the regime where the N-body results should be 
accurate. 



and likewise 



Ppow{kr) / drPpowikr)fi3{r) + Ppowikr) / di^Pwig{kr)fi3{r) 



252(27r)2 

+ Pwig(fcO / drPpow(fcr)/i3(r)+Pwig(fcr) / rfrPwig(fcr)/i.3(r) 



(27) 



where f22ir,x) and /i3(r) are short hand for the fully 
expressed kernels in Egs.lMlandl^ For the terms where 
fpow(fc) appears twice - a.k.a. the powerlaw-powerlaw 
interactions - this result can be looked up in |55| or com- 
puted using their approach. But since those results are 
often cutoff dependent and/or in poor agreement with 
simulations, we can potentially replace the powerlaw- 
powerlaw interactions and the linear theory powerlaw 
with a fitting function from pure powerlaw simulations, 
while still treating the remaining terms in P22{k) and 
Pisik) without any approximation. In Fig. 1111 this ap- 
proach is dubbed "SPT-H" while "SPT" refers to treat- 
ing the powerlaw-powerlaw interactions as in [ssj . We 
discuss the remaining interaction terms in the next two 
sections. 



E. Powerlaw- Wiggle Interactions 



Eqs. [25] and [57] contain three terms that include both 
-Ppow(fc) and Pwig{k). Since these terms include dimen- 
sionless factors of (j'o/^bao)"^^, whereas in the remaining 
"wiggle-wiggle" interaction terms there appear factors of 
(?'o/?'bao)^*""*'^\ at fixed ro/rbao these powerlaw-wiggle 
interaction terms will generally give larger corrections to 
Pbik) than the "wiggle- wiggle" interactions, which are 
discussed in the next section. The powerlaw-wiggle terms 
were evaluated numerically to obtain the SPT and SPT+ 
results in Fig. 1111 The Pi3{k) powerlaw-wiggle intcrac- 
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FIG. 12: Results for the powerlaw-wiggie interactions (cyan) 
and the wiggle-wiggle interactions (orange) for the three dif- 
ferent powerlaw setups. Expectations from the diffusion 
model coupled with an analytic approximation for Pic{k) in 
the To /r bao <C 1 limit are shown for comparison in red. Posi- 
tive corrections are shown in solid lines, negative corrections 
are shown in dashed lines. Note that for clarity the n = —1.5 
plot is shown at the first output; the n — —0.5 and —1 plots 
are shown at the second output, which allows easier visual 
comparison with the linear theory power spectrum. 
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Pwigikr) / drPpow{kr)fi3{r) 



Ppow{kr) / drP^ig{kr)fi3{r) 



(28) 

For the second term in Eq. [^Hl since Pwig(fc) is expo- 
nentially damped at high k and Pwigik) — >■ constant for 
fc — >■ 0, the result is cutoff independent. By using an 
approximation to the Pis^k) kernel one can obtain a re- 
markably accurate approximate solution for this expres- 
sion, which will be explained in the section on "wiggle- 
wiggle" interactions where this integral also appears. 

The integral in the first term in Eq. [28] also appears 
in the calculations of [s^ for a variety of powerlaws. In 
this case IR divergences might be expected to be prob- 
lematic, but, as explained by [s^l, for steep powerlaws 
the IR divergence cancels with a corresponding term in 
P22{k) (in the present context the powerlaw-wiggle term 
in Eq. 1^5)) yielding finite results for n > —3. Unfortu- 
nately, there are still UV divergences for the n = —0.5 
and —1.0 cases. We integrate up to fcmax/fcbao ~ 160 in 
the results presented here. 

The last powerlaw-wiggle interaction term, as just 
mentioned, is the second term in Eq. 1261 This term has 
a factor of two in front of it because a symmetry in the 
P22{k) kernel implies that if Ppow(fc) and Pwig(A:) are in- 
terchanged the result of the integral remains the same. 
We use this property to cross check the numerical inte- 
gration of this term. Although we were unable to find 
an approximate analytic solution for this term, we note 
that the dx integral can be computed analytically using 
the approximation Pwig(fc) sin(fcrbao)/fc and with a 
substitution of variables. This approximation is valid for 
k ^ cTj^aoi still a relatively wide and interesting range of 
k. 

The results for the powerlaw-wiggle interactions are 
shown in cyan lines alongside the linear theory spectrum 
(gray solid lines) in Fig. [121 Also shown (in orange) are 
the wiggle-wiggle interactions discussed in the next sec- 
tion. A negative correction in this plot is shown with 
dashed lines, while positive corrections are shown with 
solid lines. Qualitatively, the results indicate that the 
powerlaw-wiggle interactions are approximately out of 
phase with linear theory and push and pull the wig- 
gles in the right places to dampen out the BAO fea- 
ture. A possible exception to this is the low k correction 
for n = —1.5, but, in fact, the wide positive correction 
around fc/fcbao ^ 0.5 seems to explain the extra power 
seen on those scales in the simulation results (Fig. [T0| . 
which was not captured by the phenomenological model 
in Eq. ^ 

For a more quantitative comparison to the powerlaw- 
wiggle results. Fig. [T2| shows a model inspired by the 
diffusion behavior seen in the correlation function. If we 
suppose that the bump broadens out as in Eq. [TS] and 
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place this ansatz for cr^aoC'^o) in the phenomenologieal 
model in Eq. 1221 then in the limit where tq /^bao is small 
we expect the wiggles to evolve as 



Pwig(fc,ro)K 



Tl + 3 



sin(fcrba 



o)/2 sin(fcrbao) 
k 



2^2 { y^'^^-fc^^?^/2 Sin(fc?^bao) 

''bao / k 



(29) 

Notice that since the linear theory wiggles grow in ampli- 
tude as the linear growth function squared (i.e. Tq'^'^ ~ 
Aa"^ in Eq. [5|), the extra factor of rg^'^ in Eq. [29] makes 
this correction grow as the linear growth function to the 
fourth power. This is the same dependence on the growth 
fimction as in SPT. We plot this expectation from the 
diffusion model - essentially — fc^ times the linear the- 
ory wiggles - alongside the powerlaw-wiggle results in 
Fig. 1121 There are no free parameters to this compari- 
son; Kn takes the same value as in § IIII Ci which gave a 
good fit to the correlation function results. 

For the n = —1 and —1.5 cases the agreement with 
the diffusion model is quite good except for the caveat 
already mentioned with n = —1.5 for /c/Zcbao ~ 0.5. For 
the n = —0.5 case, the shape of -Pi3+33,pow-wig(fc) agrees 
well with the diffusion model but the amplitude is about 
a factor of four larger. Fig. [12] suggests that the prob- 
lem lies in the SPT-I- prediction, which predicts too much 
damping of the BAO feature. (Increasing the high-A: cut- 
off would predict more damping.) 

Comparing the diffusion model, which oscillates as 

— sin(A:rbao) in Fig. 1121 to the powerlaw-wiggle inter- 
actions also reveals a slight phase difference between 
^'i3-l-22,pow-wig(fc) and the diffusion model expectations. 
This is most easily visible for n = —1.5 in Fig. [12] which 
seems to oscillate as — sin(fcrbao+'rf5) where ip sa 0.2, while 
this phase is closer to « 0.1 for n ~ —1 and is consistent 
with zero for n = —0.5. This result implies that, in addi- 
tion to damping the BAO feature, the powerlaw-wiggle 
interactions provide a shift, since a Taylor expansion of 
sin(A:rbao/Q;shift) yields 

sin (fcrbao/ "shift) ~ sin(/crbao) - ("shift - 1) cos(A:rbao) 

(30) 

and, without any approximation, 

- sin(fcrbao - (/?) = - cos (/? sin(fcrbao) - s'mip cos(A:rbao)- 

(31) 

The last term on the right in Eq. [21] should provide the 
"push" to move the BAO feature to smaller scales, since 
sin (p > for the </?- values that match our results. 



F. Wiggle- Wiggle Interactions 

Though suppressed by a factor of (ro/rbao)"'*''^ rela- 
tive to the powerlaw-wiggle interactions, in Fig. [T^] the 
wiggle-wiggle interactions arc not completely negligible 




p , (k) 

wig^ ' 

Poo . . (k) 
P,, • • (k) 
P fk) 

13+22,wig-wig^ ' 
0.5 1 1.5 



2 

k/k, 



FIG. 13: Highlighting the wiggle-wiggle interactions and 
showing, individually, Pi3,wig-wig(fc) (red) and P22,wig-wig(fc) 
(blue) which destructively interfere to produce the final result, 
^3+22, wig-wig (fc) (cyan). Also shown for comparison are the 
linear-theory wiggles (Eq. 1101) in black and a cosine function 
with a similar envelope and amplitude as P13+22, wig-wig (fc). 
Note that the j/-axis is normalized to be dimensionless and 
independent of powerlaw and epoch; see text for more details. 



(at least for the n = —1.5 case), and by eye they ap- 
pear about a half-period out of phase with the linear 
theory wiggles, just the kind of feature that gives rise 
to a shift of the BAO scale. We discuss these calcula- 
tions in this section, with the convenience that because 
the functional form of P^ig{k) is independent of n, the 
wiggle-wiggle interactions are also independent of n apart 
from the (ro/rbao)^'"^"^'' term out front. Since Pv,ig{k) ^• 
constant at low k and Pwig(^) decays rapidly to zero at 
high k, the integrals should be cutoff independent. 

The task, then, is to evaluate the two remaining 
terms in Eqs. [26] and [27] We treat both terms nu- 
merically, but, fortuitously, a remarkably accurate solu- 
tion can be obtained for Pi3,wig-wig(^)- Using Pwig(^) ^ 
exp(— A:^(T^g^^/2) sin(A:rbao)/fc, and by approximating the 



Pi3(fc) kernel with /i3(r) 
488/5 one can show that 



-(352/5) exp(-29r711) 



drP^ig{kr)fi3{r) 



(32) 
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which is accurate to better than 8% for all k and bet- 
ter than 1% for k/k^^^o ^ 0.6. The minus signs in this 
result imply that, when multiplied by Pwig{k) to obtain 
^"13, wig-wig (^) as in Eq. [57] the result will osciUate like 
- sin(fcrbao)- 

In Fig. 1131 which shows the results for numerical inte- 
gration of the wiggle-wiggle interactions, the j/-axis has 
been normalized to be a dimensionless quantity that is 
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independent of the power law and epoch of interest, i.e. 

ri+3 



-^wig(^) 



^bao 



p 



bao 



^bao 



2(n+3) 



Clearly there is a great deal of destructive interference 
between Pi3,wig-wig(^) and P22,wig-wig(fc) in Fig. [I3l The 
sum of these terms, Pi3+22,wig-wig(fc)j which is of course 
much lower in amplitude than either P13, wig-wig (fc) or 

^22. ,wig-wig 

(fc), seems to oscillate at about a half-period 
out of phase with Pwig(^) as mentioned earlier. To 
highlight this we overplot with a green-dashed line a 
function proportional to — cos(fcrbao), which qualita- 
tively follows the oscillations in -P13+22. wig-wig (fc) rather 
well. Since the -P22, wig-wig (fc) term seems to oscillate as 
sin(fcrbao — v) where Lp is small and positive, when 
added to i^i3,wig-wig(fc), which oscillates as — sin(fcrbao) 
and with a similar envelope, these waves interfere as 

-sin(fcrbao) + sin(fcrbao " <<5) 

= sin(fcrbao)(-l + cos(p) - sin(pcos(fcrbao) 
f« - sin(p cos(fcrbao)- 

(33) 

The green-dashed line, more specifically, shows this 
— cos(fcrbao) term multiplied by the analytically-derived 
envelope for P13, wig-wig (A) (i-e. Eq. [32] with appropri- 
ate constants and factors of k and including a factor of 
exp(— A:^(Tbao/2) from Ppow(fc)) and divided by a factor 
of four (i.e. sin(/5 « 1/4) to approximately match the 
amplitude of -pL3+22,wig-wig(fc)- This model is only ap- 
proximate - for example, there seems to be some weak 
fc-dependence of the phase ip in P22.wig-wig(fc) - but, qual- 
itatively, something like this phenomenological descrip- 
tion must be going on. 

This raises the question of whether, in SPT, the shift in 
the BAO scale comes primarily from the phase lag in the 
powerlaw- wiggle interactions or from P13+22. wig-wig (fc)- 
The answer, at least for n = —1.5 where the BAO 
scale moves significantly, is that the shift is similar in 
magnitude from both terms, and that both "push" the 
BAO scale in the same direction. Qualitatively, the same 
can be said for the n = —1 case, but the phase lag 
in the powerlaw-wiggle interactions is smaller and the 
(''o/''bao)'"'*'^'-suppressed amplitude of wiggle- wiggle in- 
teractions is smaller still, so much less of a shift is ex- 
pected. And in the n = —0.5 case there does not seem to 
be a phase lag in the powerlaw-wiggle interactions, while 
the wiggle- wiggle interactions are even more attenuated. 



G. PT Results in Real Space 

Returning to the "coupling strength" RGPT scheme, 
which is closely related to SPT, we show the results from 
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FIG. 14: The results from fourier transforming the power 
spectrum predictions of "coupling strength" RGPT (Fig. Ilip 
into correlation functions (dot dashed lines), compared with 
the results from simulations (points). The fourier transform 
was performed with a small amount of damping in order to 
suppress noise and the influence of the power spectrum for 
k ^ fcNL, well beyond the regime where P{k) predictions are 
expected to be reliable. 



integrating the P{k) predictions from this scheme into 
two-point correlation functions in Fig. [TH Note that 
some of the outputs for the n = —0.5 case are omitted for 
clarity. At each output we apply a minimal damping to 
the power spectra to suppress noise and the influence of 
P{k) for k 3> fcNL in the final result. Some PT schemes 
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naturally include exponential damping in the predicted 
PQh{k) [e.g. [5^, which is advantageous for computing 
^(r) from PT. "Coupling strength" RGPT (and SPT) 
do not naturally include these factors, so the results for 
^(r) may not be as clean-looking as other schemes, even 
though the P{k) predictions may be quite reasonable. 
In SPT, for example, the P{k) predictions for k > /cnl 
with our setups are often large and inaccurate or predict 
P{k) < at some k. Therefore we do not show ^(r) pre- 
dictions from SPT, which offer little insight in judging 
the accuracy of the scheme or in confirming the picture 
of how the BAO feature evolves as sketched out in the 
previous two sections. 

With that disclaimer, the "coupling strength" RGPT 
predictions do a good job of rendering the evolution of 
the BAO feature in configuration space (Fig. fT4|) . In all 
cases the broadening and attenuation of the bump are 
qualitatively accounted for, including the n = —0.5 case 
that was problematic in SPT; the success of the scheme 
in this case may even help explain why the area of the 
bump is not as precisely conserved as in the other setups 
(Fig. [5]). And although it is difficult to see the trend (in- 
ferring ^(r) from P{k) over a finite k range, as described 
above, causes oscillations even when P{k) is predicted 
perfectly), in the n = —1.5 case the scheme does seem to 
accurately predict the shift in the BAO peak. With the 
close correspondence between "coupling strength" RGPT 
and SPT, broadly speaking we interpret the success of 
"coupling strength" RGPT in Figs. [11] and [H and the 
typically sensible results for SPT discussed in the previ- 
ous two sections to imply that perturbation theory can 
accurately capture the non-linear evolution of the BAO 
feature with our class of initial conditions. 



VI. DISCUSSION AND COMPARISON WITH 
ACDM 



TABLE II: ACDM outputs 



ro/'^bao 


erg 


0.003 


0.25 


0.019 


0.5 


0.040 


0.75 


0.062 


1.0 


0.106 


1.5 


0.218 


3.0 



ACDM well past z = and beyond the freeze out limit 
{as ^ 1.3). As in Fig. [Tl the overall amplitude of the 
BAO feature at fixed ro/rbao is more similar to the 
n ~ —0.5 case than to the cases with more large scale 
power. The models for the non-linear shift from [2^. 
shown with vertical dotted lines in the center and right 
panels of Fig. [TJ predict shifts of 3 — 4 % when extrapo- 
lated to the final output^. The center panel also shows 
the smooth £,n-w{r) correlation function, computed from 
a fourier transform of Pnw (k) from (60j , and in the right 
panel ^nwir) is subtracted from the simulation data. In 
the center panel the combination of strong damping of 
the BAO feature and noise in the ^(r) measurement make 
any shift non-discernible. In the right panel the result of 
subtracting out ^nw(?') docs visually resemble an atten- 
uating gaussian (much more than S,(r)/^nwi'f), which is 
not shown), but it is unclear whether the apparent drift 
of the BAO peak towards smaller scales, especially by the 
last output, is truly from the non-linear shift or whether 
the effect is simply from the changing broadband shape of 
^(r). A plot of (^(r) — ^pow{r)) / {z) versus r from any 
of our fiducial simulations would show a similar trend. 



A. ACDM-like Simulations 

Having described and explained the non-linear evolu- 
tion of the BAO-feature with our powerlaw setup in some 
detail, it is worth discussing the relevance of these results 
to the canonical ACDM cosmology. We approach this 
task first by simply assessing the resemblance of our re- 
sults to ACDM. To aid in this comparison we performed 
a set of four simulations with an initial ACDM spectrum 
{n,n = 0.226, Oa = 0.774) as in Fig. [2] but evolved with 
flm = 1, ^^A = SO that tTg and ro/rbao in this case 
can avoid the freeze out limit and reach values compara- 
ble to the powerlaw setup. The ACDM-like simulations 
presented here were performed with essentially identical 
parameters as the earlier fiducial simulations in terms of 
box size, force resolution and number of particles. We 
show the primary ^ (r) results in Fig. [151 the ro / rbao val- 
ues for each output is shown in Table [HI 

Fig. [15] is fairly unremarkable except that it shows 
the non-linear evolution of the correlation function in 



B. Perturbation Theory and Modeling 

In§ IV Bl we showed that a phenomenological approach 
matched the results from our fiducial simulations rather 
well. Eq. [22] bears a close resemblance to the damped- 
exp oncntial models often used in the literature [e.g. 
[22l^3j. and we emphasize our conclusion that the broad- 
ening (damping) of the bump (wiggles) depends on the 
pairwise dispersion, Spaij., rather than the rms displace- 
ment. S^, which is sensitive to bulk motions. In Fig. I16[ 
we compare Epjjjj./S^ on a wide range of scales for a 
ACDM spectrum (Fig. Although we expect the two 
formulae to converge to the same result as r — >■ 00, it 



The prediction depends on whether one assumes their ashift — 
1 oc D{z)^ formula, as expected from SPT, or instead uses their 
empirical fit where aghift — 1 tx D(2:)^ '^^. Fig. 1151 shows the 
predictions of the D{z)'^ model. The empirical model is similar. 
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FIG. 15: Left panel: The matter two-point correlation function results from four simulations using a canonical ACDM linear- 
theory power spectrum but evolving the initial conditions using Sim ~ 1, ~ 0. The correlation function is shown at 
different epochs (points with error bars; solid when ^(r) > 0, circles when ^(r) < 0) with the linear theory correlation function 
overplotted (solid colored lines). The vertical dotted line shows the initial mean interparticle spacing. Center panel: The 
correlation function near the BAO scale. Vertical dotted lines show the expected shift from [SJ] colored according to epoch. 
Also shown is the smooth ^nv,{r) (black dot-dashed line) derived by fourier transforming Pnw(fc) from [gO]. Right panel: The 
result of subtracting ^nw(r) from the ^(r) measurements. 




pairwise displacement 

---S2 = l/37r2/„°°/'i,(9)d9 



10" 



FIG. 16: The rms pairwise displacement (Eq. 1111 solid) at 
different scales using ACDM initial conditions divided by the 
commonly used rms displacement formula (Eq. O dashed), 
which includes the contribution from bulk motions. Both 
quantities scale as the linear growth function squared, so the 
result shown is independent of epoch. 



is nevertheless surprising that I]p3^jj,(rbao) differs by less 
than 2% from the displacement. In the literature 
some groups use Eq. [2]to predict the damping, while for 
others is a free parameter that is fit to simulations 
[e.g. m. In our view, like that of [H, it is E2^i^(rbao) 
that matters physically, and the success of models based 
on Eq. [2] is a lucky coincidence that holds in ACDM-like 
models but can fail, by an infinite factor, for powerlaw 
models. 

Another widely-used phenomenological approach as- 
sumes a model for PneC^) motivated by Renormalized 
Perturbation Theory [RPT;[6l|. In these models the non- 
linear shift comes directly from including P22{k) in the 
phenomenological form, or, in real space, from modeling 
the shift with the closely- related ^mc(?') ansatz and cali- 



brating the amplitude of this term to N-body results (e.g. 
psl l29l l45j). Using our setup and a natural value for the 
amplitude of this term, in § IIII Dl we showed that this 
approach adequately captures the shift in real space for 
the first output of the n = —1.5 case (when a% = 0.5). 
By the second output (corresponding to as = 1), how- 
ever, it fails, and although not rigorously justified by the 
derivation of the term, we argue that the formula would 
more accurately predict the shift if the broadening of the 
bump could be incorporated into ^nic{f)- This may have 
been previously unnoticed because the shift in ACDM 
when cTg ~ 1 is smaller than the shift in the n = —1.5 
case, and the amplitude of the bump, i.e., ^(rbao), is sig- 
nificantly smaller in ACDM than in the n = —1.5 setup. 

Finally, the success of the "coupling strength" RGPT 
method [l[ in matching our simulation results, both in 
fourier space and in real space, may certainly be infor- 
mative to ongoing efforts to model the BAO evolution 
with ab initio predictions from PT. [l^] show that this 
scheme also does a reasonable job in predicting the non- 
linear power spectra of ACDM and cCDM cosmologies. 
Except for SPT [s^] we ignored other PT schemes, but 
in principle the predictions from many other PT schemes 
could be compared to our simulation results and useful 
insights gained from the kind of comparisons presented 
in § IV CI This would no doubt be useful for BAO studies, 
and, more broadly, find that the largest deficiency of 
the halo model is in capturing the transition from the 
1-halo to 2-halo term, precisely the scales where the per- 
turbation theory predictions arc most important. 



VII. SUMMARY 

Motivated by the importance of accurate modeling of 
the BAO feature in large scale structure for interpret- 
ing the results of future dark energy experiments, we 
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have used N-body simulations to investigate the evolu- 
tion of a BAO-like feature in a simpler, alternative set- 
ting, where it modulates an underlying powerlaw initial 
power spectrum in an ilm = 1 universe. Specifically, our 
initial conditions have the correlation function defined 
by Eq. [l] with a gaussian multiplicative bump centered 
at scale rbao and the amplitude ^bump and width CTbao 
chosen in approximate agreement with ACDM expecta- 
tions. The corresponding initial power spectrum follows 
Eq. [TU] to an excellent approximation. For given values 
of ^bumpi fbaoi and the powerlaw spectral index n, non- 
linear matter clustering statistics (including the correla- 
tion function and power spectrum) should depend only 
on the ratio ro/rbao; where rg is the correlation length 
defined by ^(rg) = 1. We evolve our simulations to values 
of ro/rbao much higher than traditional ACDM models, 
with final outputs corresponding to ag ~ 4.0 {n = —1.5), 
6.0 (n = —1), and 12.0 (n = —0.5) if we define a phys- 
ical scale by setting rbao = 100/i^^Mpc. Our standard 
simulations have box side Lbox/fbao = 20 and 512'^ parti- 
cles. We use our simulations to develop physical intuition 
for BAO evolution and to test analytic descriptions in a 
regime far from that where they have been tested previ- 
ously. In this respect, the spirit of our exercise is similar 
to the cCDM investigation of ^ and [U.^ 

Consistent with ACDM studies, we find that the 
strongest effect of non-linear evolution on the BAO fea- 
ture in ^(r) is to fiatten and broaden the bump, with 
^bump decreasing and crbao increasing. To a good ap- 
proximation, failing only at late times in the n — —0.5 
model, the area of the gaussian bump, proportional to 
^bump X Cbao, remains constant, which suggests that pairs 
are "diffusing" out of the shell corresponding to the ini- 
tial BAO feature (see the physical description of f2^). 
The evolution of the bump width is well described by 
a model in which the non-linear crbao is the quadrature 
sum of the initial width and a length proportional to 
Spain the rms relative displacement (computed from lin- 
ear theory) of pairs separated by r = r-bao- The constant 
of proportionality varies with n, but the same constant 
that describes our standard n = — 1 model also describes 
the faster evolution of an n = —1 model with a "skinny" 
initial bump, supporting the validity of the diffusion in- 
terpretation. For n = —1.5 (where the relevant integral 
converges without a small scale UV cutoff) the diffusion 
constant computed ab initio describes the bump evolu- 
tion accurately. We emphasize that it is Epair rather than 
the rms absolute displacement S that is relevant to ana- 
lytic descriptions of our models. The latter quantity has 
an infrared divergence for n < — 1. but this divergence 
corresponds to bulk translations induced by very large 
scale modes, which cannot affect the BAO peak itself. We 



Another notable study is [63l ] who investigated the non-linear 
evolution of a ACDM spectru m p lus a fourier-space spike on 
scales relevant to BAO. [g] and [25l | have also discussed toy mod- 
els for BAO but without investigating non-linear effects. 



think that the appearance of E rather than Spair in many 
analytic models of BAO evolution is at best an approxi- 
mation restricted to CDM-like models with a turnover in 
P(fc); by coincidence, E « Epair(7'bao) for ACDM. 

The location of the BAO peak, defined by the scale 
''peak of a gaussian fit to the non-linear ^(r) divided by 
the linear theory powerlaw, stays constant within the sta- 
tistical precision of our measurements for the n = —0.5 
and n = — 1 models, even when these are evolved to a 
highly non-linear stage where the bump amplitude has 
dropped by a factor of ~ 4—10 from its initial value. For 
n = —1.5, on the other hand, the peak location shifts 
to smaller r, an effect that is already noticeable at the 
first output (rg/rbao = 0.024, equivalent to erg = 0.5) 
and that grows to a 30% drop by rg/rbao = 0.386 (equiv- 
alent to as = 4.0). The analytic models of [H] and 
accurately predict that shifts should be much larger for 
n = —1.5 than for n = —0.5 and n = —1, and the [25| 
model accurately describes the evolution of the peak lo- 
cation for n = —1.5. However, both models predict non- 
linear shifts in the n = —0.5 and n = —1 cases that are 
inconsistent with our simulation results at late times. 

We carried out a number of additional numerical 
tests varying either numerical parameters or the phys- 
ical model. Our fiducial simulations have ibox/^bao — 20 
and an initial mean interparticle spacing smaller than 

— 1/3 

''bao by a factor of rbao/'^p =25. We found con- 
sistent results in simulations with L^ox/fhao = 10 and 

— 1/3 

Thao/np ~ 50, indicating that a box size ten times the 
BAO scale is acceptable. We found marginal discrepan- 
cies for 256^ simulations with rbao/'^p = 12.5. Success 
of the box size test and other internal consistency tests is 
achieved only because we include the integral constraint 
corrections described in Appendix [B] which make a no- 
ticeable difference for n = —1 and an important differ- 
ence for n = —1.5. In other tests, we show that BAO 
evolution is nearly identical in an = 1 model and 
a model with ri,„ = 0.3, JIa = 0.7 (and the same ini- 
tial conditions) provided they arc evaluated at the same 
value of rg/rbao (or, equivalently, the same value of the 
linear growth function). 

For more thorough tests of analytic models, we turned 
to a fourier space description using the non-linear matter 
power spectrum. A "phenomenological" model in which 
we combine numerical results for the non-linear power 
spectrum of a pure powerlaw model (Appendix [X] and 
references therein) with our gaussian fits to the evolu- 
tion of the BAO bump in ^(r) gives a remarkably ac- 
curate description of the full non-linear outputs of the 
n = —1 and n = —1.5 models. This model assumes no 
shift of the ^(r) peak location for n ~ —0.5 and n = — 1 
and rpoak/j^bao = 1 - 1.08(rg/rbao)^'^ for n = -1.5. The 
success of this model suggests that the BAO bump has 
little effect on the non-linear evolution of the underlying 
"smooth" power spectrum. At least for rg/rbao < 0.2, 
we expect that this model is a more accurate descrip- 
tion than our numerical P{k) measurements themselves, 
since it draws on self-similar scaling results from pure 
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powerlaw spectra that have wider dynamic range than 
our simulations. 

We compared our results to predictions of two ab ini- 
tio analytic approaches, "standard" 1-loop perturbation 
theory (SPT; e.g. iHsll) and the "coupling strength" 
RGPT scheme of [l| . This scheme provides a quite accu- 
rate description of the low-A: evolution in all cases, includ- 
ing n — —1.5 where the peak location shifts significantly, 
and it produces good but not perfect agreement with the 
evolution of the ^(r) bump in configuration space. For 
SPT, we break up the terms into distinct "interactions" 
between the powerlaw and "wiggle" components of the 
linear power spectrum, both to obtain physical insight 
and to allow us to define a more accurate "SPT+" scheme 
that uses numerical results for pure powerlaw evolution 
and perturbation theory to describe the interaction terms 
that involve the "wiggle" spectrum. SPT alone gives 
a reasonable description of the early P{k) outputs for 
n = —1.5, but on the whole "coupling strength" RGPT 
is substantially more accurate and has a wider range of 
validity. 

The high statistical precision achievable with future 
BAO surveys — with ^ 0.2% cosmic variance distance 
scale errors for z > 1 and redshift bins Az « 0.2 [64 1 — 
puts stringent demands on theoretical models. Exploit- 
ing the power of these surveys will require large numeri- 
cal simulations supplemented by the physical insight and 
modeling flexibility afforded by analytic methods. The 
simulation results presented here offer valuable "stress 
tests" of numerical and analytic approaches in regimes 
beyond those where they are usually applied, and they 



allow isolation of distinct physical effects. Two natural 
directions that we plan to explore in future work are the 
clustering of biased tracers — in particular the massive 
halos expected to host luminous galaxies — and the im- 
pact of redshift-space distortions on BAO measurement 
from galaxy clustering. We will also investigate the im- 
pact of the initial conditions algorithms, comparing the 
scheme advocated by [65j for simulation ensembles to the 
traditional scheme of mean density boxes used here. The 
combination of future BAO surveys and improved the- 
oretical models will lead, ultimately, to new insights on 
the energy and matter contents of the cosmos. 
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Appendix A: Results from Pure Powerlaw Simulations 



Having performed pure powerlaw simulations for the sake of better understanding the non-linear power spectra 
of our fiducial simulations, we give fitting functions for the n = -0.5,-1 and —1.5 powerlaws using 512'^-particle 
Gadget2 simulations, which were set up similarly to the fiducial simulations as outlined in §|TT1 The interested reader 
can consult the excellent paper by (37| to find fitting functions for other powerlaws. 

Fig. 1 171 shows our primary powerlaw results compared against other fitting functions in the literature, either specific 
to each powerlaw or universal fitting functions designed to match a variety of p owerlaws and cosmologies. Our 
simulations do not extend to impressively large values of fc/^NL compared to [37[, in part because of how long we 
chose to evolve the simulations and in part because we chose, conservatively, to only show /c-values up to one fourth 
the particle nyquist wavenumber, i.e., half the rule of thumb recommended by |52| . However, we run a number of 
realizations of each powerlaw (six realizations for n = —0.5, four for n = —1, and ten for n = —1.5), which is a few 
to many more than in previous studies. As a result the error bars in Fig. 1171 which show the measured errors on the 
mean from all realizations in each case, can be quite constraining. 

Since there is always a concern that the numerical results will be invalidated when the clustering power on the 
scale of the box becomes large, following the convention of [13 we show the value of a/a* = (fcB/fcNL)''"''''^-'^^ for each 
output in all three panels as an indicator for how close the non-linear scale has come to the scale of the box. As stated 
previously, even our simulations with the most large scale clustering (n ~ —1.5) fall comfortably below the threshold 
where the loss of clustering power from beyond the box scale might be a concern. More quantitatively, in Fig. 1171 
a/a* is typically ^ 1, and in the n ~ —1.5 case the last output only reaches a/a* = 0.18. Importantly, the later 
outputs seem to show the self-similar scaling required by the scale-free nature of the initial conditions, the results 
falling along the same curve when plotted against fc/fcNL and divided by A|(A:). For the earlier outputs this scahng 
seems not to have set in yet in some cases, a fact revealed by the self-similar test. Therefore we define the non-linear 
fitting functions as much as possible to the later outputs which arc least affected by the clustering signature of the 
initial grid. 
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FIG. 17: Results from pure powerlaw simulations (colored points with error bars) compared to various fitting functions in the 
literature (black dotted, dash-dotted or dashed lines). Also shown is our improved fit in thick black lines (Eqs. [XT] & [A2| and, 
for n = —1.5, the analytically derived SPT prediction from Appendix B of [S^] is shown with a thick black dashed line. For the 
Smith et al. (2003) functions in the n = —1 and n = —1.5 panels, we use modified formulas that correct typographical errors 
in the original paper (see eg. IA3l and following text). 



We present our non-linear fitting functions as a generalization of the functional form in [37 

A^L(fc) = A2(A:)/„(A:/A:nl), 



(Al) 



f l + Ax + Bx" + Ex' 
■'"^''^ ~ 1 + Cxf + Dx' 



(A2) 



Our n = —1 results primarily drive the necessity of making this generalization. The n = —1.5 results seem 
reasonably well represented with the Widrow formula, so we set £' = £' = £ = 5 = in that case, while in the 
n = —0.5 case we set £ = e = 0, which still allows sufficient degrees of freedom to adequately describe the simulation 
results. The fitting formula above should be accurate to k //cnl '-^ 5 for ri = —0.5, k/kj^n^ ^ 100 for n = — 1 (larger 
because the fit was forced to closely match the results of [33] at high fc), and fc/ZsNL ^ 6.5 for n = —1.5. 
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TABLE III: Best-fit Parameters for tlie Non-linear Fitting Function 



n 


A 


B 


C D 


E 


Q 


P 


7 


5 


e 


-0.5 


-0.1309 


0.1131 


0.1296 -0.02472 


0.0 


8.599 


2.066 


8.714 


0.4565 


0.0 


-1 


-0.4722 


0.3542 


0.04449 -0.2020 - 


0.08956 


1.358 


1.447 


1.911 


0.3963 


0.2564 


-1.5 


-0.0792 


0.1704 


0.008748 0.0 


0.0 


1.225 


2.672 


2.1306 


0.0 


0.0 






1. 


Specific Comments 


on n = 


-0.5, 


-1, & 


-1.5 







To our knowledge the most recent work to explicitly show the non- linear evolution of a pure n = —0.5 spectrum 
from simulations is l35i. whose universal fitting function we plot alongside our results in Fig. 1171 The universal fitting 
fimctions of [6^ and j36j are capable of making predictions for n = —0.5, but their fitting functions were trained only 
on ?7 = and n = — 1 simulations to set the scaling in this regime. With this caveat the remarkable agreement of the 
prediction of [66j and our n = —0.5 results seems somewhat fortuitous and the disagreement with the [36j prediction 
seems not so surprising. 

The discrepancy between our n = — 1 simulation results and the Widrow fitting function is entirely explainable by 
the sparseness of their measurements in the quasi-linear regime and is not indicative of any kind of problem with 
either their or our simulations. At larger /c/Zcnl both results overlap nicely, and we define our fitting function to 
closely match theirs for /c/Zcnl 3. Also plotted alongside the n = — 1 simulation results is a fitting function specific 
to n = — 1 from Appendix B of |36|. There is a typo in their fitting formula (their Eq. Bl), which should read 



fsdsiy) = V 



1 + y/« + {v/hf + [y/cY 



(l + (zj/d)("-'3)7)l/7 



(A3) 



(R. Smith private communication). The corrected formula for = — 1 is shown in the middle panel of Fig. 1171 and 
although at fc ~ /cnl it deviates strongly from either fitting function, at lower k it matches our results reasonably well. 

Appendix B of [s^ also includes a set of constants tuned specifically for their n = —1.5 results, but there is a typo 
in their table in the reported value of a. From quantitative comparison to their Fig. 11 (especially at large k/kjsii^), 
the correct value seems to be a sa 7, rather than a = 0.707 as rep orted. We show this result alongside our other 
results for n — —1.5 in the right panel of Fig. [iTl Since both [6Q| and [s^ include n = —1.5 simulations in their 
universal fits, we also show the predictions of their fitting functions. Finally, we plot the expectations from SPT for 
n = —1.5 in the limit that the UV cutoff goes to infinity, as in Appendix B of |55| . In our fiducial simulations this is 
the formula used to predict the evolution of the powerlaw in the SPT model shown in Fig. [TT] The SPT+ models in 
Fig, [m instead use the non-linear fitting functions just described to model the evolution of the n = —0.5 and n = — 1 
powerlaws. 

Appendix B: Integral-Constraint Corrections to the Measured Matter Autocorrelation Function 

When estimating the correlation function in our simulations, we divide the average number of neighbors found 
around particles in the separation range r — > ?■ + rfr by the number expected for an unclustered distribution of number 
density N/V: 

i{r) = (Nn^^i- ^ - + dr)) ^ ^B^^ 
' Airr^dr x N/V ' ^ ' 

where V is the simulation box volume, is the total number of particles, and we use ^(r) to distinguish this estimated 
correlation function from the true correlation function ^truo (f) of the underlying cosmological model. This procedure 
is subject to a well known "integral constraint" bias (described by, e.g., [s^ §47]), which arises because the simulation 
volume itself is forced to have the cosmological mean density. The fact that the total number of particle pairs in the 
box is N{N — l)/2 « (1/2)7V^ imposes the requirement 

fls=ibox/1.61 



C(r) « / inr^dr^ir) = , (B2) 

Jo 

where we have approximated the integral over the box volume as the integral over a sphere of volume (47r/3)i?| = 
V = i^ox- large volume ACDM simulations, the bias in ^(r) is usually a small effect because the true correlation 
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function goes rapidly to zero, then becomes negative at large r, making equation (jB2j) easy to satisfy. However, for 
powerlaw models with negative n, the slow decay of the correlation function makes the integral constraint bias more 
important. 

We account for the integral constraint by assuming that it produces a scale-independent additive bias, so that the 
mean value of ^(r) averaged over an ensemble of simulations would be 

C>) =etmc(r)+€bias . (B3) 

For our powerlaw models, Eq. (|B2p then implies 

f-Rs 

47rr2dr [CtrucW+ebias] = (B4) 







and thus 



where we have used the linear theory ^L(r) = (r/ro)"^"'*''^' for ^truc(?')- More elegantly, this bias is simply the volume- 
averaged correlation function, ^bias = —^L(i?s), which agrees with the conclusions of [63, §6.4.2], who derived this 
term using the sophisticated error analysis in [68|. 

In all our figures we plot the corrected correlation function 

^ i{r) + URs) ■ (B6) 

At large r. the fractional correction is 

^(r)-^"(r) 3 /1.61r 



n+3 

(B7) 



Since r is always less than Lbox/1-61, this correction is fractionally larger for more negative n and, at fixed ri, the effect 
is most important for r approaching the box scale as previously mentioned. In practice, we find that the integral 
constraint makes little quantitative difference to the appearance of, e.g.. Figs. [3] & [8] for n = —0.5, a noticeable 
difference for n = —1, and an important difference for n = —1.5. In particular, the box size convergence tests in 
Figure [5] succeed for n — —1.5 only because we include the integral constraint correction. 
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